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Abstract 

A  procedure   for   the  direct  numerical   integration 
of  the   steady-state,    elastic   scattering   neutron  trans- 
port equation  is   presented.      Conditions    for   the   ex- 
istence,   uniqueness   and  non-negativeness    of  the   nu- 
merical solution  are    obtained.      Under  these   conditions 
it   is   shown  that   the    scheme   is    stable   and   the   numerical 
solution  converges,    with  refinement   of   the    mesh,    to   the 
solution  of  the    transport   equation. 
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Sectioh  1«      Introduction 

II         ■    ■        <  ■ ■  .■■■—■■■->-—  ^  H  I    I   ■■  ■■■MM 

A  procedure  for  the  direct  numerical  integration 
of  the  one-space-dimensional  transport  equation  has 
be6n  derived  and  tested*   The  success  of  thefee  comp- 
tations  has  made  a  general  description  of  the  method 
desirable.   Thus  the  particular  case  of  a  mult i*- region, 
infinite  plane  shield  is  considered  here  as  an  illustra- 
tion of  the  general  technique. 

The  method  is  practical  for  steady  state  problems 
in  which:  1)  diffusion  theory  is  not  ^plicable,  2) 
slowing  down  by  elastic  scattering  is  important, 
3)  finite  geometry,  internal  boundaries,  or  varying 
composition  are  present.   Anisotropic  or  inelastic 
effects,  which  permit  the  integration  of  the  birth  in- 
tegrals in  section  2,  may  also  be  included.   If  either 
the  spatial,  angular  or  energy  dependence  may  be  elimin- 
ated the  method  becomes  practical  for  various  reactor 
calculations  (i.e,  for  a  particular  Fourier  component 
of  the  flux  in  a  slab  reactor) , 

The  finite  difference  equations  and  an  explicit 
solution  procedure  are  given  in  sections  i|.  and  5«  Here 
sufficient  conditions  on  the  mesh  spacing  are  obtained 
to  insure  that  the  numerical  equations  have  a  xinique, 
positive  solution.   For  practical  computations  it  is 
believed  that  these  conditions  may  be  relaxed. 
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In  sections  6  and  7  we  show  that  the  numerical 
scheme  is  stable  and  convergent*   This  means  that 
"small"  errors  introduced  at  the  boundaries  influence 
t)^   solution  only  by  a  "small"  amount;  and  by  taking 
a  sufficiently  fine  mesh  the  numerical  solution  can 
be  made  as  "close"  as  we  please  to  the  true  transport 
solution,   I'he  conditions  under  which  this  is  true  are  the 
same  as  those  which  are  sufficient  for  the  solubility  of 
the  equations.  Although  the  convergence  character- 
istics of  any  method  are  of  central  importance  in 
determining  its  practicality,  this  property  is  gen- 
erally taken  for  granted,  without  any  proof.   This 
matter  is  treated  in  section  7;  the  method  of  proof 
also  furnishes  a  bound  on  the  error  which,  although 
pessimistic,  may  in  sorre  cases  be  useful. 

Section  2,   The  Transport  Equation 

The  steady  state  transport  equation  for  neutrons 
with  elastic  scattering  and  isotropic  source  (i,e, 
linearized  Boltzmann)  is 


(2.0)  [7*-a+  (5^(r,u)]^(r,n,u)  =  Hb^^^l  *(r,i^,u)  J  + 

+  — S(r,u)  ; 
where 
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2.1   b^^hv'(r,/L,u)]  s  \         du'<5f(r,u') 
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r^t'  ,u-). 

A  discussion  of  this  equation  is  contained  in  [l]  and 
the  notation  used  is  as  follows: 

Dl   r  =  position  vector  of  neutrons  in  the  shield  or 

reactor; 
D2  /I  =  unit  vector  in  the  direction  of  neutron 

velocity; 
D5   u  =  In^  "^^^^E,   lethargy  variable  defined  in 

terms  of  the  maximum  source  energy; 
D4   i/'(r,/l,u)  =  neutron  flux  (=  vN(r/l,u),  where  v  is 
the  neutron  speed  and  N  is  the  neutron  density 
in  (r^,u)  space); 
D5  6    (rju)  =  total  cross-section; 

c 

D6   fi>(r,u)  =  scattering  cross-section  of  the  i-th 
scatterer; 

D7   A.  =   atomic  weight  of  the  i-th  scatteror; 

^       A  +1 
D8  1±  Z   21j^  a  .1  >  maximum  lethargy  degradation  for 

neutrons  on  collision  with  nuclei  of  the  i-th 


scatteror; 


u-u'        u'-u 


D9   a^(u-u')  =1  (A^+1)  e  2   -(A^-l)e  2 

of  the  scattering  angle  in  the  lab  system  for 
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elastic  scattering  from  the  i-th  scatterer 
with  lethargy  degradation  u-u' ; 

DIO    F^^^(n./ii.u,u'  )  = 
(A.+l)^ 


U^A^ 


p.  (a.(u-u»),u«  )e^  "^6[A.Jl«-a,(u-u»  )j. 


probability  that  a  neutron  scatter  from  (A',u')  to 
(A,u)  on  collision  with  the  i-th  scatterer  (for  low 
energies,  u»»0,  and  frequently  for  light  scatterers 
it  is  usually  assumed  that  p.  (a.,u»)  =  /2); 
Dll    S(r, u)  =  flux  of  source  neutrons. 

The  sum  over   i   in  (2,0)  is  assumed  to  include  all 
elements  present  at  position  r.   The  source  term,  which  is 
discussed  in  the  next  section,  is  assumed  known.   This 
essentially  restricts  the  present  considerations  to  shield- 
ing problems,   A  further  restriction  is  imposed  upon  the 
geometry  and  symmetry  of  the  problem  as  follows:  define 


^2.2)      ^.  =  -L./^, 


the  cosine  of  the  angle  between  the  position  vector  of  a 
neutron  and  its  velocity  vector.   It  is  then  assumed  that 
the  flux  has  the  functional  dependence 

(2.3)        ^(r,A,u)  =  ^(r,ti,u)  . 

This  assumption  is  valid  in  at  least  two  important 
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cases:  first  infinite  plane  gfedm^trles  in  which  all 
spatial  variations  depend  upon  one  coordinate,  say  x; 
and  second,  spherically  symmetric  geometries  in  which 
all  spatial  variations  depend  upon  the  radius,  r.  In 
the  plane  case,  which  will  be  the  main  consideration 
of  this  paper,  (2.0)  becomes 

(2.4)  [ti^  +  6'^(x,u)]^(-x,ti,u)=  T:b^^h*(x,ti,u)]  + 

+  ^5^S(x,u)  . 
For  spherically  symmetric  geometries  we  have 

(2.5)  tt^j^  +  ^  ^  +  <3^(r,u)]vl;(r,ti,u)  =  Jh^^hll{r,^,u)]   + 

+  ^(r,u)  . 

In  either  of  these  cases,  or  in  general  whenever 
(2.3)  is  valid,  the  "birth  integrals"  (2,0)  may  be  simpli- 
fied by  explicitly  carrying  out  one  of  the  angular  inte- 
grations.  This  can  readily  be  done  if  we  introduce  the 
spherical  coordinates 

n  =  (1,0,0)  >i  ^ 

(2.6)A'=  (cose»,sin©'cos0»,sin©'sinj^»  )  (  ,  )  ^   _,     > 
r  (p<0,0'<2:^ 

-rV  =  (cos©,  sin^cos^,  sin©sin0)      I 

in  (2.1)  to  obtain,  using  (DIG)  and  (2,3), 
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(2.7)    b^^h^]    =i   i  du'g.  (r,u,u' )    j    d0'    (de'sine'SL 
max   4 


cosO'  - 


For  simplicity  we   have    introduced 


-a^(u-u'  )]\|/(r,  |^',u«  ). 


(D12)    g.(r,u,u')    =  —i dj(r,u')p,  (a,,u»)e''  -''     , 

^  '^  ""     2t[A^  ^  ~ 

and   from    (2.2)    and    (2.6) 


^iv-^, 


(2,8)    [I  -  cose        , 

tx'    =   cosQcos©'    +   sin9sin©' cos  (j^' -0)    • 

A   standard  theorem  on  5-functions    states 


5[g(c)]  =  n 

k 


dg(?w) 


dc 


-1 


5(4-€i,)      , 


where    the    sum  is    over  all   the    zeroes    (assumed   simple) 
of  g(4);    and    it    is    simple    to   verify   that 

2%  271  71 

r  f(cos[0'-0J)d0'    =   rf(cos0")d0"   =   2  f  f(cos0'"  )d0'" 
^o  ^  ^o 

Using  these  resu'^.ts  and  (2,8)  in  (2.7)  yields 

u  71 

(2.9)  b^^^[\l/]  =/du«g^(r,u,uM/d0\Kr,^^u') 


(,U-q. 


where 


(2,10)  r^^   =   §^(ti,u,u«,0)  =  |aa^(u-u«)   + 

+  [  (l-lJL^)(l-a^(u-ut))]   2cos0  . 

For  purposes  of  the  numerical  treatment  it  is  useful 
to  introduce  the  integral  operator 
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(2.11)  h^^^[\l;(r,ii,u);u«]  =  fd^^Kr,  S,,u«  ) 

Where  S^  is  given  in  (2,10).   The  birth  integrals, 
(2,9),  can  now  be  written  as 

u 

(2.12)  b^^^[\l/(r,^i,u)]  =  f  du'g^(r,u,u')h^^^[4,(r,ii,u))a'j  . 


max 


[u-q^ 


Section  3.      Plane   Syirmietrical   Shield 

Consider  an   infinite  plane    isotropic   source   of 
neutrons,    infinites imally  thin,    perpendicular  to  the 
X-axis   at   the   origin.      Then  the   source   term,    (Dll), 
is   given  by 

(3.0)  S(x,u)    =  5(x)f(u)    , 

;^?here  f(u)  is,  say,  the  fission  energy  spectrum  of 
neutrons  with  f(u)  =  0  for  u<0  and  f(u)>0  for  u>0 
(see(D3)).   Let  this  source  be  surrounded  by  a  finite 
symmetric  shield  in  the  range  -a<x<a  ,  with  composition 
a  function  of  x  only.   Outside  the  shield  we  assume 
vacui;im  of  infinite  extent.   This  is  obviously  a  case  of 
infinite  plane  geometry  in  which  condition  (2,3)  holds 
with  r  replaced  by  x.   The  appropriate  form  of  the 
transport  equation  is  then  given  by  (2,1^),  with  the 

source  (3,0),  and  the  birth  integrals  (2.9)  in  which 

T 
r  is  to  be  replaced  by  x,  V/e  allow  6     and  g«  to  be 
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piecewise   continuous    fiinctions   of  x;    thus    the    shield 
may  be   composed  of   layers   of   different      composition  and 
variable   density. 

At   the  vacuum-shield   interfaces    there   can  be  no 
neutrons  entering    the   shield,    thus 

\l;(a,ix,u)    =   0   , -l<n,<0  , 
^^•■'•^      \l/(-a,ti,u)    =   0,l>ti>0   . 

From  the   symmetry  of   the   problem  it   is    clear  that   the 
flux  miist    satisfy 

(3.2)  \|/(x,n,u)    =  \l/(-x, -tj,,u)    . 

With  this  result  it  is  possible  to  eliminate  the  source 
from  {2,l\.)    and  replace  the  above  problem  by  one  in  the 
range  0<x<a,   Thus  we  integrate  (2,ij.)  with  respect  to 
X  from  -e  to  e,  apply  the  symmetry  condition  (3.2), 
and  pass  to  the  limit  e  — >  0,   Since  all  spatial 
singularities  in  the  cross-sections  and  fliix  are  assumed 
integrable  this  yields 

(3.3)  v|'(0,ti,u)  =  \l/(0,-|x,u)  +  i  f(u). 

The    problem  has   now  been  reduced   to   solving 

(3.i+)      ttiT^  +  <3'^(x,u)]\l/(x,^l,u)    =  j;b^^^[^(x,ii,u)]    ; 
subject   to   the    boundary  conditions,    from    (3.1)    and    (3«3)» 
a)\l/(a,ti.,u)    =   0,    -l<[x<0  , 

(3.5)  "  1 

b)\i/(0,p.,u)    =  \!/(0,-ij.,u)    +  ^  f(u),    0<|x<l   ; 
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in  the  domain  (see  fig, II) 

(3.6)   Dy?  [0<x<a,  -l<ia<l,  0<u<U] 

At  the  maximum  source  energy,  u  =  0,  the  birth 
integrals  (2.9)  vanish  and  the  flux  of  highest  energy 
neutrons  must  satisfy,  from  (3, if), 

(3c7)  [\i-^   +  6^(x,u)mx,[x,0)    =  0  , 

and  the  boundary  conditions  (3*5) •   The  solution  is 

easily  shown  to  be 

^0  ,-l<li<0 

(3,8)\l/(x,ti,0)  =^  ,  0<x<a   , 

X  "~ 

f(0) 


i^exp[::ijd^(C,  0)d^],  0<^<1 


The    singularities   at  n.  =   0  in    (3.5)b)    and    (3.8)    are  due 
to   the    infinite   extent    of  the    source   plane   which  conse- 
quently has   an  infinite   fliix  of  neutrons    streaming   along 
its   surface.      These   singularities    are   untenable   for   the 
proposed  numerical   procedure    and  may   be   eliminated  by 
introducing  the    "slightly  anisotropic"   source 

/5(x)   Mf(u)    ,    0<|ix|<tx. 
(3.9)      S(x,ti,u)   =S 

^6(x)f(u)     ,  K;,<!ti|<l 
Here  cos"  [j,,,.  is  an  arbitrarily  small  angle.   This  source 
differs  from  the  isotropic  source  (3.0)  only  in  a  small 
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angular  region  where   the    isotropic   distribution  is 
replaced  by  a   cosine    distribution.      The  boundary 
conditions   for    this   source   become 

\|/(a,  iJ,,  u)    =   0  ,    -l<\i<0 

(3.10)  (i  f(u)    ,    0<^<^., 

(if(u)      ,    lx,.<tJ.<l 
The   highest  energy   flux,    whiclr  satisfies    (3.7)    and    (3. 10), 

is   given  by 

\|;(x,p,,0)    =   0  ,-l<^t<0 

il/(x,ix,0)    =   exp[-i/  6Ui,0)da)Y'-^ 

r-f(O)    ,ii,,Si<l^ 

The   effect   of    the    "slightly  anisotropic"   source   is   shown 
in  fig. I  v/here   the   highest   energy  flux   in  the    source   plane 
is   plotted   as   a    function   of  \i  for  both    (3.8)    and    (3,11), 
The   spatial  dependence   of   \};(x,n,,0)    may  be    inferred  by 
considering   exponentially  damped  waves   emanating   from  the 
appropriate   \{/(0,  tj,,0)    curve    o£    fig. I,      The   rate   of   damping 
must   increase    as    [x  — >  0  through  positive   values. 

In  summary  we   seek  the    solution    (assuming   existence 
and  uniqueness)    of   equation    (3»i|)    subject   to    the   boundary 
conditions    (3. 10)    in  the   domain  D^  of    (3,6).      The    flux 
at  highest   energy  in   this    solution   is   given  by    (3,11), 
Prom  physical   considerations    the    solution  should  be   non- 
negauive   and    continuous    in  x.      The   domain  D„   is    shown   in 

-13- 


!'      ■■■% 


'    '  '     (,'. 


'-  -\ 


'f  J  ■> 


o.- 


-<)■.,■■     .  *  f  ■  *■■  .•  ^  • 


/'■;     /' 


'■''/■.■•. 


■J- 


'<.'} 


'    ■■'    '...I 


fig, II  where    the    appropriate  boundary   conditions    and 
highest    energy  solution  are    indicated. 

Section  li.      Finite   Difference   Approximation 

In  this    section  the  problem  formulated  at   the 
end  of  the   previous   section  will   be    approximated  by 
a   set   of   finite    difference   equations.      Thus    on   the    domain 
D^j  of    (3 •6)   we    introduce   the   mesh/i^r  defined  by  the   points: 

^)^r+l  "  V-^^r  ^    ^o  "^  °  *    ^R  "   ^'    (0<r<R); 

(i|,0)        ^    -^  s        s  o  o^  ..^+1  «2.  *^ 

°'^n+l  ""  %  "^^  V    %  "^   °  '    ^N  "^  ^  '    (0<n<N); 
d)^j^l  =  0.   +A0y    j2^o  =   0  ,    ^j  =  71   ,     (0<j<J), 

Here  the  0.    are  values  of  the  dummy  Integration  variable 

J 

in    (2.9)    or    (2,11),      All   spacings    4  x^,    ^[i^,    A '^^  and 
^  0.   are   positive.      The    angular  spacing    (lj.,0)b)    is 
chosen  such  that   M-o   t^  ^   ^°^   ^H   ^»      This   eliminates    the 
consideration  of   a   special  set   of  difference   equations 
for   the    case   ix  =   0,      The  x-spacing,    (i|,0)a),    is    to  be 
such  that   all    interfaces    (points   of   discontinuity  of  the 
cross-sections)    coincide  with  mesh  points.      Values    of 
functions  which  are   known  at   points   of  ^^j,  are   denoted  by 
a)(5^(r,n)    =  6^U^,^^)    ,    see    (D5); 

b)g^(r,n,nM    =  Si(x^,u^,u^,  )    ,    see    (D12); 
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c)f(n)  2  f(u^)  *  see  (3.0)  and  (3.9)/  , 
(k.l) 

d)S^(s,n,n',  j)  =  ^i^l^s*  V^"^j^»  ®®®  (2.10), 

T 
If  X  is  a  point  of  discontinuity  of  6  (or  g.)  then  (L|..l)a) 

or  b)  is  defined  as  follows: 

'^^(^r-0,u^)  If  0^(r-l,n)  Vpp^.^s  in  the 
(i|..2)6  (r,  n)  =  )  (difference  equation 

■"  %  rp       (with  6T(r,n); 

6-'(x^+0,u^)  if  d^(r+l,n)y 

(and.  similarly  for  g.).  This  lack  of  s ingle- value dness 
of  the  cross-sections  on  /!>.  allows  the  limit  of  the 
numerical  solution  to  be  continuous  and  have  discontinuous 
derivatives  in  x.   The  meaning  of  (I;. 2)  is  made  clearer 
after  (Ij..li|). 

With  some  aim  at  simplicity  in  the  present  treat- 
ment all  integrals  are  approximated  by  the  trapezoidal 
rule  arxi  only  linear  interpolations  are  used.  Before 
applying  these  techniques  it  is  useful  to  consider 
fig, III;  a  typical  (^L,u)-plane  at  some  fixed  value  of  x, 
say  X  =  X  ,  in  which  the  domain  of  integration  for  the 
1-th  birth  integral,  (2,10)-(2.12),  is  indicated  by  the 
shaded  region.   The  intersections  of  horizontal  and 
vertical  solid  lines  are  points  of  the  grid  (ij.,0)b),c) 
for  the  fixed  x  .   The  dots  represent  locations  of 
f .  (s,n, n' ,  j)»  which  function  as  the  |x-argument  in 
\|f(x,^,u)  (see  (2.11)),  for  some  fixed  (s,n)  and  all 
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(n»,j)  i^L  0<j<J>  n.<n'<n.   Here,  as  indicated  in  the 
figure,  n.  is  defined  as  the  integer  such  that 

u  <Li,  -q.<u  ^,,  for  u  -q.>0  , 
1         1 

n.  =0        ,  for  u^-qi<0  • 

Let  V/^        be   the   numerical   approximation   to   the 
r,  s 

flux  \|/(x^, p.  ,u   )    on    A^,      Then  the   integral   operator    (2.11) 
is    approximated   on  the  mesh  by  the   difference   operator 

^^^^^^^,s''^'^-^^^^^^^V^'s'%^'''n'^    *   "^^^"^^ 
a)H^^Uv/5^3;n.]   =iZ^^f[,3      1.5^(s,n,nSJ)]l/5;3      -. 

-f   [5i(s,n,nt,j).^.3^3w5]3^^l|,    n^<n' <n; 

(l+.U)  ^  ^ 

b)H^^hw5^3;n]    =  itw5^3      . 

Here   s  .    is   the    integer  such  that-::- 

ik»S)      l^.,    <^i(s,n,n«,  j)<ti3   ^^      , 

J  J 

and  we   have    introduced 

,    J   =   0 

(i|.6)      e^   =  3^i2^j.i  +^i2fj    ,    l<j<J-l 

The   definition    (l|,tj.)b)    follows    from  the   fact    that 

a.  (u  -u   )    =   1   and  hence,    S.(s,n,n,  j)    =  [x      for   all    J, 
inn  1  s 

There   is    obviously  a   dependence   of  s.   on  i,n  and  n»    which, 
for   simplicity  of   notation,    has  been"^ dropped.      This    is    true 
also  of  the   sets   L   ,    and  R^,    about   to  be    introduced. 
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For  practical   numerical   procedures    and  for  the   stability 
analysis    it    is   liseful   to  write    (l4.,l4.)a)    as 

(k.7)      H^^Uw^     :n']   =21    ^,(s,n,n',s')V/^\,    ,    n.<n«<n. 
r,  s  s'=0   ■■■  ' 

Prom   (i|, 0)d)    and    (ij..i|) -(i|.,7 )   we  must   have 

(I|..8)    u^(s,n,n»,s«)   =  MZ     J^^^^^s  +i-S^(s»n,nt ,  j)  ]   + 

where  the  disjoint  sets  L  ,  and  R  ,  are  defined  by 

L   ,    =  <  all    j   such  that:    u,„      =  M-   , 

(k.9)  ^  ^  ,    0<s'<S2 

Rg,    =  >  all   J  such  that:    tig  ^^^  ="  l^s' 

Clearly,    from  the   definitions    {k*5)    and    (L|.,9),    we  have 
in  the  set   theoretic   sense 

which  aids    in   confirming    the    identity  of   the   sums    (Ii..l4.)a) 
and    {k»7)»      If  some    of  the    sets     ^Lg,+Rg,^      are   empty, 
which  may  frequently  be   the   case,    the   corresponding 
^  .  (s,n,n',s' )    vanish   (see   fig. Ill), 

The   birth  integrals    (2.12)    are   approximated  on 
^jj  by  the   difference   operators   B^^^wJ^  ^3  ^b^^h\|<(x„,iig,u^)  ], 
which  are   defined   as 
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ik.lO)    B^^hv.^^J    =1     ^5.(n')g^(r,n,n«)H(^hv;5^3;n»],n>0. 


Here   n.    is   defined   in    {k»3)>    H^^^V/^  „;n^  ]    is    defined   in 
([;,i|)-(ij..6),    g^(r,n,n«)    is    defined  in   (if,l)b),    and  we 
have    introduced 

An-1  ,n'    =  n  ; 

I^%«-1   "^^^n'  ,n>n'>n^   +   1   ; 

/  1 

1  "n  ^1     n<      2 

Uu     [l-( i)]  ,nt    =  n,      . 

1 

In  order  to  construct   a  solution  of  the   equations  being 

formulated   and  to   analyze   their  stability  it    is    useful 

to  write    (i|,10)    in  the   form 

(k.l2)    B^^^lwJ^J    =^u^^3_  |g^(r,n,n)w5^3   +  B^twJ^s^'    ^>°   '. 
where  we  have   used    (i4.,l4.)b)    and  the   definition 

(t^.l3)    B^[W5  J    =i|:     6    (n«)gi(r,n,n«)H<^)[wJ     .-n'l,n>0. 

It  has  been  assumed  here  that  ^u  _^-i<q.  •   If  this  is  not 
the  case  4  ^v,-!  ^^   Ci^lS)  should  be  replaced  by  the 
0.(n.+l)  of  ([j.,11).  However  it  is  felt  that  the  present 
procedure  should  not  be  used  (at  least  in  its  current 
form)  in  this  case  and  thus  (L|.,12)  will  be  retained 
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throughout. 

The    difference   equations,    which  are   to   approximate 
the   transport  equation    i3»k)    o>^4vr  are   chosen  to  be 

"t-l't'  1  ,.,      „  ,.,  (^0<r<R-l 

i  (    l<n<N 

The   difference   quotient  is    centered  at  x     i/,  =  x     +  — ~ 

and    the    other  terms   are   taken   as   averages   about   this   point, 

T 
If   there    is    a  discontinuity   In  6     or  g.    at  some   x  =  x     it 

T 
is   now   clear  that   the    value    to  be   used  for   6    (r, n)    or 

g.(r,n,n')    is   that  belonging   to  the   medivm  containing 

the   point   about  which  the   appropriate   equation   (ij.,ll4.) 

is   centered    (see    (i|..2)    and  following).      To   these   equations 

we   adjoin  the    approximation  to  the  highest  energy  flux 

(3.11) 

w°  s  =  0  »o<s<S3_; 


f(o) 


(i^.15)  ..1..     ,_  _  .    lUT'^i^  ^  >• 


r 


.   =    (exp[-i     i:    ^l6^(r'+l,0)-f6^(rt,0)|])]  i^Si<t^s^M 


Ii2l,aii  s   3: 
^^s 


for  0<r<R  , 
The  boundary  conditions  (3,10)  are  approximated  by 
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Here  we  have  defined  s_  as  the  largest  integer  <S_  such 
that  M-o^-M-o* 


Section  5«   Solution  of  the  Difference  Equations 

In  this  section  an  explicit  procedure  for  solving 

the  difference  equations  (l+.ll^. )-([;.  16)  is  presented. 

This  provides  a  constructive  proof  of  the  uniqueness  of 

the  numerical  solution  which  is  required  of  any  practical 

scheme  and  is  necessary  for  the  stability  analysis, 

Sufriclent  conditions  on   the  mesh  spacing  are  presented 

which  insure  that  the  solution  procedure  is  possible  and 

that  the  numerical  solution  is  non-negative. 

We  first  note  that  equations  (l4.,li;)- (1^,16),  with 

(4.7)  and  (1+,10),  form  a  system  of  (N+1)  (R+1)  (S2+I) 

inhomogenous  linear  equations  in  as  many  unknowns,  V7*,   , 

r,  s 

The  solution  of  this  system  is  obtained  by  an  induction 
on  n.   To  start  the  induction  we  note  that  the  highest 


energy  approximation,  W°  ,  is  given  in  (l\.,lS)    for  all 
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r  and  s.   Then  assuming  given  the  t/J   for  all  n»<n-l 

r,  s  — 

and  0<r<R,  0<s<Sp  we  write  ([[.ll;)  as 
[tXg+^ix^p(r+l,n)]w5^^^3-[tx^.Ax^p(r,n)]w5^g=-^(B^^[w5^^^J 

Here  we  have  used  (i|.12)  and  (ii, I3),  and  introduced 

(5.1)  2p(r,n)  =  c5^(r,n)  -  Ziu^.;^  I  2Igi(^»n,n)  . 

The    inductive  hypothesis    insures,    by   (k»7)  and(l4.,I3), 

that    the    right    side   of    (5.0)    is    known.      For  s<S, ,    l"/?        =   0 

by  the  boundary  conditions    (ij.,l6)a).      Then   starting  with 

r  =  R-1,    and  continuing   recursively  down  to  r  =  0,    equations 

(5.0)    yield   all   w"   ^    for   0<r<R-l   and   s<S, ,      The    values    of 

W^        so    computed  yield,    by    (ij.,l6)b),    all    the   Vr^        for 

s>S, ,     Another  recursion,    from  r  =   0  to   r  =  R-1,    in    (5.0) 

yields   the    remaining  W  for  l<r<R,    s>S- •      This   completes 

the    induction. 

The    recursive   solution  of    (5.0)    is   possible   if 

a)  -tig   +^x^p(r,n)    j^  0     ,      3^^' 

(5.2)  (  0<r^-l 

b)  [i^  +>^x^p(r+l,n)   7^0,      s>S^ 

and  the  induction  holds  if  these  conditions  are  satisfied 
for  all  n>l.   With  the  much  stronger  conditions 
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iB.3)^  /  ^  ^0<r<R-l 

it  is  easily  proven  that  w"  >0  on  the  entire  mesh  <^,t. 

r,  s—  N 

This  proof  relies  on  the  fact  that  w  (s,n,  n» ,  J)>0,  and 
f(n)>0  for  all  i,n,n',s  and  j.  The  first  of  conditions 
(5«3)a)  and  b)  are  satisfied  provided 


iS,k)     i3(r,n)>0  ,   0<r<R  ,  n>l  . 

These   conditions   imply    (5»2)    and  hence   the    solvability 
of  the    difference   equations.      To   examine   condition    (S*k-) 
in  more  detail  we   use   the  definitions    (5«1),    (D12)    and 
(D9)    and  recall  from  the   definition  of  total   cross-section 
that 

(5.5)        o'^(r,n)    =XIt<(r»n)+6?(r,n)]. 

i      ^  ^ 

This    yields 

(5»   P(r,n)   =^^65(r.,n)^[l-A^„.,-^  p^(l.n)  ]uf  (r,  n)^  . 

where   6f(r, n)    is    the    absorption  cross-section  of  the    1-th 
scatterer  at    (x^,u^)    and  p^(l,n)   =  Pi(cLi (^^'^n^' ^n^  * 
The    lethargy  spacing  may  now   be   taken  such  that 

6^(r,n)  kA. 

(S.l)  4u„   T<[l-»-4 ] o^ *^or  all   i,r  and  n>0; 

^-■^-       6^(r,n)      (Aj+l)'=^p^(l,n) 
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in  which  case  (Sik)    is  satisfied*   When  this  condition 
is  satisfied  the  second  of  (5»3)a)  ^^^   h)  are  Satisfied 
provided  the  spatial  and  angular  spacing  are  chosen  such 
that 

(5»o)  ^  X  <  ,  for  all  r,  s  and  n>0, 

inax[p(r,n),p(r+l,n)] 

Conditions  (5»7)  and  (5.8)  are  then  sufficient  for  the 

existence  and  uniqueness  of  a  non-negative  numerical 

solution,  l/^  ,  on  the  mesh  -d„.   For  practical  compu- 
r,  3  vi  ^ 

tations  these  restrictions  are  no  doubt  too  severe. 

Section  6.«   Consistency  of  Difference  Equations 

In  the  present  section  we  outline  a  proof  of  the 
consistency  of  the  difference  equations;  that  is,  that 
the  difference  equations  approach  the  transport  equation 
as  the  mesh  is  refined.  A  trivial  observation  of  con- 
sistency in  the  present  case  is  obscured  by  the  complexity 
of  the  equations.   The  present  analysis  also  yields  a 
bound  on  the  truncation  errors,  which  is  required  in  the 
convergence  proof. 

Let  \|((x,  ^i,  u)  be  any  function  satisfying  (3,10) 
and  (3,11)  with  differentiability  properties  in  D  which 
are  to  be  specified  as  required.   Then  we  introduce  the 
linear  transport  operator  '?'[  ]  by 
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where  b^^h\|(]  is  defined  by  (2.11)  and  (2.12).   The 
difference  equations  are  represented  by  the  linear  differ- 
ence operator, 

(6.1)  i I ♦  (x^^^,. u„)  i^J'l^Y^,  s^  ^^(r*l,n)*;^^_^  * 

where   \|(^      ==  \k(x  ,  ti,  ,u   ),the  mesh  is     zl^  defined  by   (i+.O), 

and   the   B^^h\|/]    are   defined   in    (k'k)    and    (ii.lO).      The 

T 
values   of   ^    (r,n)    (and  g^(r,  n,n'))    at   points   x  =  x^  of 

discontinuity   are   defined  by    (i|.2).      We   define   the 

truncation,    T^  „,    on    A.r  by 
r,  s  JM 

(6.2)  (r-t)[yl^(Xj.^^,ti3,u^)]   =  Tj^g        , 

for  the   given   function  \|>.      The    difference    operator,    t, 

is   said   to   be   consistent  with   the   transport   operator,    p^   , 

for    the  given  \|/   if 

(6.3)  T^  g — >  0  as  inax(4x^,  ^tig,  ^  V^^j^ — ^  ^  °^    "^  N  ' 
The   proof  proceeds    from  an  evaluation  of   the    truncation 
by  standaird  methods. 

From  the   usual   Taylor  expansions  we   obtain 
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(6.1^) 

P        P 

2     (^yx 
Here  we   have   used   the    obvious   notation,    ^'j.^.ta  g  z. 
\l/(x  +  ,  LL  ,u   );    and  the    assumptions   that  <±i^     and 

*^— '^ — r^     are   continuous    functions    of  x   in  each  interval 
<^x2 

[x    ,,-x    ],      The   contributions   to   T^  „    from  the    birth 
r+l     r  r,  s 

integral   approximations   are   obtained   from  the   identity 

The    superscript    i  has   been  temporarilly  dropped  for 
simplicity.      As    in    (6,Lt-)b)   we   now  have 


(6.6)  blt?^y2^^1-|[b[t5^,_J+blt5_3i]  = 

=  — -^— 2.  b[*(x^+e4x^,^,,u„)]  , 

<-  ffX 

with  proper  continuity  assumptions  on  «— w  •  The 

<)x 
remaining  terms  in  (6.5)  are,  by  (2.12)  and  (i+.lO), 
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(6.7)  b[^^^3J.B[t5^3]=2:     ^    du'rg(x^,u^.u.)h[^5^^;u.). 

-  •|(g(r,n,n»+l)H[4/J^g,n'+l]+g(r,n,n')H[>|/^^g,n'l)  ?, 

In  order  to  avoid  a  special  case  of  no  particular  diffi- 
culty we  have  here  assumed  u  -q.  =  u   ,  The  above 
integrand  can  be  rewritten  as 
f  -  ?=  6(x^,u^,u')h[i|»5^g;u'].  ■^(g(r,n,n«+l)h[^^^g;n«+lj  + 


(6,8) 


+  g(r,n,n«)h[il/^^;n«j) 


+  |g(r,n,n»+l)(hl^)/^^g;n«+l]- 


-H[^^^3;n«+1]) 


+  ^(r,n,n«)(h[y>^^g;n«]- 
-H[4'S  .;n»]) 


An  expansion  similar  to    that   of    (6.l|)b)    (but  with  respect 
to   u'    in  the    first   two   terms   on  the   right  of    (6,8)    yields 


u_,+ 


(6 


nVX. 


•9^    J   du'5g(x^,Uj^,u»)h[rl/5^g;n'J-  |(g(r,n,n»+l)h[il/J^3;n'+l]  + 
^n»       C 


+g)r,n,nMh[^l;5^3,n«])U 

-   ^  ^  V  )^  c)g(  W^*  )h[^g^  3;u'  3 
1^  c)^' 


U'ssU 


n' 


+er(/au^, )  . 


The  remaining  terms  in  (6,8)  are  evaluated  from  the 
definitions  (2,11)  and  ([|.,li.)-(i^,6)  which  yield 
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:::.:V./-V.- 


'r.^-  -,::^ 


'^-v^-( 


r  » 


— -^  \ 


^^.r.r-i^^). 


:'..S.  .K- 


,;■■  ^'•':?'i^.».,.  •,. 


■'-'•-    ~,-r,. 


■■  ;  '1 


y  J  ?n^  *.-.): 


•<•■■  •  •• 


•^^.•'■:i       ,  :.. 


;vJ..;> 


'  '  •  v* 


;  L 


.■t\.  •■,•. 


-■'.//?.■.' 


J-1  0-^^0y 
(6.10)  M^l/^  3;n'].H[rI/J  g;n']=2;  fd0'UfU^,^(^,^,u^,u^,jS),n^,) 


J-1 


l^s.+l-^^'^s'W^j) 


j=0   ^  ^V-s  J 

ffi(|Xs,u^,u^,,jZfj)-^ig  J-1    , 

^  — ^;: ^^^^j-^r  ^^o'^^jf  v^^^s' vv^i^H..)- 


^l^. 


^'^J-Hl 


'j  +  1 


A\^, 


J+1 


■r,sj^^+l- 


The  usual  Taylor  expansions  are  now  used  above  to  obtain 


(6.11)    h[vl/^      :n»]-H[4/2  .;n«  ]  = 


r,  s 


J-1 


r,  s' 


0=0. 


0=0. 


-o'Ca;^^)^ 


J+i 


J-1 


^^t^s  )^  ^^s,+l-^^^^s*W'^J)  ^^^^s'W^j)-^^s.'^^?,s•. 


k 


i-.[-Ii 


^l^c 


A  tie 


c)M- 


ZllJ-c 


■J+1 


^ti< 


'J+1 


c)l^ 
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•       if  '  <  ;  : 


U.-> 


7-:  U, 


W 


■-!f'  ^V 


V>.N      


'  :>(•;■ 


■n    ■!■■;-;■ 


^    r  .I.:  '\  ■-  I  ,-■-    ••,;■■    •  .■>.'' 


l^^> 


r.  V  ■■■ 


-k.. 


»  . .    :  ^  ■•'  "   '  4- 


■.;l 


C-.:. 


!  V 


Here  we  have  introduced  \|/^  —  =\|/(x„,  ti^. +©4^.  ,u  )  for  some 

r,Sj-       r     Sj  Sj.     n 

0<0<1    «        Using    (6,[i-)-(6,ll)    in    (6,2)   we  may  obtain  an 

expression  for  T^      •      Since   this    result   is   rather   cumber- 

some   we   present    instead  a  bound   on  the   truncation  which 

is    obtained  from  the   complete  expression.      Thus   we   let 


( 6 ,  12 )  A  x=max Ax  ,A  u=max An  ,  Z^ii,=max4 tx  _ , 4  ^max60 .  j 
~r       ^       "n       "       ~s        ^       ""j        J 


and  then 


(6.13) 


T 


r.s 


<(4x)^a+(4ti)^P+(^u)Y+(^i2f)5+»(^u^+4|2f^) 


Here   a,    p,    y   ^^^  ^    ^^^    independent    of   the   mesh  spacing 
and  depend  only  upon   the   cross-sections    and   particular 
fiinction  \|/   and  their  derivatives   in  D^,     Explicitly 
^e  have   defined 


,2^T 


a  =  -^j-imax  ^^--*f-»-max  "  -    J^-  +ymax 


P  =  f  |ZI  "^i  "^^^  gj^-max 


^^b^^U.] 


c)x' 


ai 


^ 


^u 


(6.1U) 


^KfK.. 


max 


Dy     ()tl 


«)g^h^^^[^l/;u] 


c)u 


5  =     ^iZI  ^.i   ^2LX  g .  »max 


c)H'(x,S^,u)  <jg^(u,u,n',jZ^) 


«^2. 


^i^ 
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{;i^)- 


■'A:f.r: 


•    !  ,i.';Li^ 
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■  ■    .  (. 


.    /   I 
..JJ\ 


.  •     .•    ^ 


'\  ,  t  V    :;    '■:'•■  .^'i- 


J  •*  ..» 


;ti^v:..;; 


It   is    clear  from   (6.13)    and    (6.1if)   that    (6.3)    is 
satisfied  provided  the    indicated  partial  derivatives 
remain  finite. 

Section  7.      Stability  and  Convergence   of  Numerical   Solution 

Let   ^{x,\i,\i)   be   the   unique   solution  of    (3,1^.), 
satisfying    (3.10)    and    (3.11),    in  D^.      Let  W^  „  be   the 
linique   numerical  solution  of    (t|..ll4.),    satisfying    (i^,l5) 
and    (i|..l6),    on  A-^,      The    error   in  the   numerical   solution  is 
defined,    on    ^^.r,    as 


<7.0)        «?,s   =<,s    -   ^^W^n)- 

If  we  define  W  }/    in  any  reasonable  way  (say,  = 

|[wj  g+wJJ^^  g])  we  have,  from  i^.k),    ik-^k) ,    (6.0)  and  (6.1) 

?-[\l;(x,n,u)]  =  0   ,    on  Dy  ; 

(7a) 

From  the  linearity  of  these  operators  and  the  definition 
(7«0)  we  obtain,  using  (7.1) 

*f^?.yp,33=t[Cv,3].t[Mx,4i,^.^,u^)] 


{r-tlin^, 


<ax 

^,ix„,u„)]   • 


^r'  2  ""s'^n 
From  (6.2)  the  above  becomes 
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v\ 


1 1    '■-•        •■-*- 
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where  the    truncation  is   estimated  in   (6,13)    and   as 
outlined   in  the   previous    section  may    be  obtained  explicitly 
in  terms    of    the    solution,    t,    and  the  mesh  spacing. 

The    errors    on  the  boundaries,    x  =   0  and  x  =  a,    are 
easily  obtained  from   (7.0),    (14-.16)    aid    (3.10)    by  the  methods 
of  the   previous   section.      We   include  here   an   inequality 
which  results    from  the    complete   expressions: 


e"       =  0 
r,  s 


,    0<s<S     ,    n>0; 


(7.3) 


e"^g|    <  max    I  ej^  ^  |    +    (Am,)^^   ,    S^<s<S2  ,    n>0; 


s<S^ 

where    A[i.  is   defined   in   (6.12)    and  ^  =  -^fmax 

u>0 


^\l/(0,n,u) 


c)u 


Similarly,  for  the  error  in  the  highest  energy  solution, 
from  (7.0),  ik.lS)    and  (3.11) 


r,  s 


,       0<s<S^      ,       0<r<R; 


(7.i|) 


,0 
'r,  s 


<    (4x)t9  ,      S^<s<Sg    ,       0<r<R; 


where  Z\x   is   defined   in    (6.12)    and 


af(0) 


max 


\   -      2K:,^s^+i   0<x< 


+  ©'(ax) 


n 


To  obtain  bounds  on  the  errors,  e  *    we  consider 

r,  s. 


-30- 


•  \ 


•.-.V) 


-.:.o) 


{-)■:: 


o'y 


(7c2)    for  r!<S-c^"Tn   this  T'ang'^--i^-.:C-^~.t-T.:o   ^^i;tain,    u:>:i.n5- 
the    derinition3    (6.1},     '.k^l2)    and    (5^1;^ 

^'^°^''         °r,  S     "     7~^,"~r'77:"~*~""r-;.i..o 


:/2 (-::;-' ;^^)r^:^. 


Let   us    introduce    EOir.e   norms    of    tho  error   by 


L„  „  =  ma;c| '^.    i  .  ^.,  ,  •   rr.axie^' 


r,n  -^  .  .A  r.  s'  -   r,.i 


-  J  'J 


'    0< d <2  ^  '  -      -  '  •-   S ^  <3 <S  - ' 

(7.6: 

\  =  ^-qej  ^.1  -  max  (L   /R^   )   . 
0<r_<:^    '     0<r<R   '    "' 
0<3<S2 

Then  from  the  definitions  (i;rl3),  (14-.  11)  and  '(l|o7)  we  havc 

—  „_  n-l 

(7  ,7  )  E:(B  '  "■  h  el,^^  3  ]  >'b^'  ^ '  [  e;j_  ^^  1  )  <^uG.^i:  \,       ' 
^  '  n-^O 

where  ^u  ir.    defined   in    (6^12 )    and  we  hf.ve    introduced 

(7o8)      G     ^  2r.  na>:  X!s.  vr,n,n^' }       e 

0<n~'<n-l 
From    (7  o  6) -(7.8)    in    {7'<5)   we   obtain;,    talcing   maxima  with 
respect   to  s. 


^7<=9)  L^,n^^n^r-.:., 

.=31- 


-r  p  : 
n.         n 


v' 


with  the  definitions: 


max  I  ||ig|-4x^p(r+l,n)  | 

(7.10)  \=^^ 

"^  mini l^^l+dx  B(r,n)| 


s 
r,  s' 


p  = ^^ 

n  —  „.. 


mln|  |^g|+4x^p(r,n)| 


n-1 


Tn^^^^l^?,sl^n^^^^nS:  \.  ^\ 
^'^  n«=0 

Applying  (7.9)  recursively  yields 

B  r       l-^^n)"""* 
^7.11)   Lr,n  ^  ^^n)   %,n  "^  " "7 ^n^n' 

In  an  exactly  analogous  manner  we  obtain,  by  con- 
sidering (7.2)  for  s>S^, 

(7.12)   R    <  (\)\^^   +  TT^^'nV 
whe  re 


max 


ltig|-4x^p(r,n) 


(7.13)  7l„=^^^-^— ■ — ,P  = 


^x/^ 


n—  . 
min 


ltig|+/lx^p(r+l,n)|   "^  min|  ltig|+4x^p(r+l,n)j 

r,  s 


r,s 

The  first  term  on  the  right  in  (7.12)  indicates  the  propa- 
gation of  errors  introduced  at  the  boundary  x  =  0.   A 
similar  remark  applies,  in  (7.11),  to  the  errors  introduced 
at  the  boundary  x  =  a.   These  errors  would  actually  be 
damped  if  \x    \<1   and  \/[    |<le   The  errors  d-ue  to  truncation 
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'■':  <.. 


■"  -i. 


.,  r 


.^ .:    :■   *■! 


\     ( 


4       •      / 


;>frv;--:. 


i> . ..-  i 


and  those  introduced  at  the  mamimum  energy  are  propa- 
gated in  the  terms  containing  E  • 

From  (7,3)  we  deduce,  using  (7 .6), 


%^,n=  0 


(7.1i^) 


Now    (7.11)    and    (7.12)   become 

i-(x  )^-^ 

L-,  „<[ ^ P„JE 


r,n- 


l-\ 


(7.15) 


n 


n     n 


>R 


R 


1-(A    )^^        l-(\)'  _  2 


Prom   (7.10)    and    (7.13)   we   see   that   in  a  shield  of 

piecewise-constant  composition;   X     =A     =  X"     and 

n         n         n 

p^  =  P^  =   p^,    say,    and    (7.15)   becomes 

L 
(7.16) 


1.(1  )^-^_ 
^  t Pn^^n 


1-X 


n 


1  (X  )^'^^ 

Rr,n  ^  t-^-^i ^PJE^  -H    (XJ^(4^)24    . 


1-X 


n 


In  general   it   is    easy  to   show   that  the   above  equations 

are   valid   if  X_   and  "p     are  defined  as 
.  n  "^n 

maxjl  |Hg|-4x^p(r+l,n)|,  [  1  t^g  | -^ix^p  ( r ,  n )  | 

r  zJLlLL 


\  ~~  >  "-^ ^  >niax(X   ,^) 

'^"minf||ixJ-f4x^P(r,n)|,||iXg|+,lx^P(r+l,n)|>  " 


r,s 


(7.17) 


r,s<J 


Ax/. 


lt^J+^.P(r,n)L||tx J+Zix  6(r+l,n)|?=^^^Pn'^n^ 


-33- 


-i^ao'ia 


.    .  ..  t  >  / 


-■--  i. 


H     V.CiO> 


xi. .  V} 


v;-'',!-    o. 


r^'ia.;:  ■  ■•■;   ■,)'i'; 


1      : 
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•■;  i   i      f^  '■"' 
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Si^ 


i  1     :  r 
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Taking  the  maximum  in  (7.16)  with  respect  to  r  and 

recalling  definitions  (7.6)  yields  an  inequality  of  the 
form 

\  <  ©n^n  ■*■  ^^^^^  ^   raax[l,  (\^)^1, 

where  _ 

e  =  max  -%-[i-(r  )R-^  ,  i-cx^)^-"')  . 

^  0<r<R  1-X      ^  ^ 
n 

Without  loss  of  generality  we  may  assume  X  <1  in  which 

case  these  relations  become 

a)   M^©n^n  "^  ^'^^^^^      * 

(7.18)  _  pp 
1-(X^)2R_ 

From   (7.10)    in  the    above    inequality  we  have 

"n^^^^n^nt:     \'    "^  Vn  ^  ^^^^     * 
n'=0 

or 

n-1 

(7.19)  M^<^t.K^Z:     M^,    +  F^       , 

n'=0 
where  we  have  introduced: 


a)        K  =  max©  G 
n'<n 

(7.20) 

With  these  definitions  we  easily  obtain  from  a  recursive 


':.)r     r\ 


application  of  (7.19) 

(7..21)    M^(l+4iiK^)''"^[^i^K^MQ+F^] 


A  well  known  Inequality  states  that 

(7.22)     (l+|)''<e^ 

for  k>0  and  all  positive  integers,  n  «   Thus  if  the 
lethargy  spacing  wore  equal,  say  n^u  =  u„»  we  would  have 
by  the  above  inequality 

{l+4uK^)^-^<e  "  '^-^ 

For  arbitrary  spacing,    as   used  in  the   present  scheme, 
there   exists   for  each   value   of  n<N   a  bounded  number, 
t     jT»    such   that 


(7.23)       4u  =  -B  t^^^  ,    l<tn,N^^' 


where  we    recall  -^u  =  max  4u   .      Thus   in  general 

'"0<n<N  ^ 

(1.2k)  (l+4uK^)'^-^<e^"^-^*''^-^'^ 

and  (7.21)  Implies 

(7,.25)    M^<e  "^   ""•*•  "^--^'^[AuK^Mq  +  F^} 

From  n ,l\.)    and  (7.6)  we  have 

(7.26)     Mq<(4x)  ^l 

From  (7,.20)b),  (7.10)  and  the  estimate  (6.13)  we  have 


■  --    -.'*  v^ 


■i'-. 


\  : 


^fi.;';;':    ;•     ;'.>Xi 


i.:ix>    f;'-c;?ii.:; 


•(  v.r- 


n 


:  cr>.       j>.i 


:■.■■■{'.> 


■.\,i-Lr:-cr:r^^.j.-L^^ 


■!'> 


;■  \ 


it 


i'-.,\)     -^':-i 
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',■■,»  1  . 


:.() 
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•  •  \  ■■■!'■ 
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•     i.vr 


(7.27)  F^efcax)2a+(dti)^P+(/lu}Y+fci^)5+e'((4u)^+(^j2f)2)  ^ 


+  (^l^)^c  , 

where 

(7.28)  Q   =  max  0 

"  n<N  "^ 

Thus  the  final  form  of  (7.25)  may  be  written  as 

K  u   t 

(7.29)  M^<e  "^  ^-^  ''•^'^[(du)(^x)K^r^+©(^x)2a+(4ix)V 

+  t<3u)Y+(AiZf)6+er((/iu)2+(^^)2)j+(^^)2^j  ^ 

We  now  assvune  that  the  mesh,  ^^,    satisfies 
conditions  (5.7)  and  (5.8).   Then  the  inequalities  (5.3) 
hold  and  the  absolute  value  signs  may  be  dropped  in  the 
definitions  (7.17).   ^ince  the  denominators  in  X  and 
p  are  evaluated  for  the  same  values  of  r  and  s,  say 
r» ,  s' ,  we  have 

Pn  ^^ 

(7.30)  — S-  = 


1"\^       2^tig,+^^^,p(r»,n)-tig„+^x^„p(r",n) 

.  I ) 

where  r",  s"  are  the  integers  for  which  the  numerator  in 
X  is  a  maximum.   For  simplicity  we  have  written  the 
P(r, n)  at  r  =  r'  and  r".   These  could  be  replaced  by  the 
values  at  r  =  r»  ■+-  1  or  r  =  r"  +  1,  as  the  case  may  be, 
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with  no  change  in  the  following.   The  conditions  imposed 

on  the  mesh  insure  that  X  <1  and  hence  for  such  mesh 

n 

spacing  the    denominator  in   (7 .30)    never  vanishes. 

The  methods   used   to   derive    (7.214.)    yield,    when  applied 
to   the   definition   (7.17) 

(7.31)  l-(5:„)^^<l-«xpf.[-^;A.,   Rp(r",n)   + 

■^  T--^^>  j^p(r",n)][2aj 

Here   r* ,    r",    s'    and  s"   are   defined   in  the   previous  paragraph 
and   as    in    (7.23)   we  have    Inti^duced  the   »cundod  pceltlve 
constants      ?^,    p  such  that: 

(7.32)  4V=t^r.,R  ' 

Condition   (5.7)    insures   that    the   p(r»,n)    are   >   0  and   thus 
t^e  exponential   in   (7.31)    is   bounded  for  all  possible 
Hg,    .     From   (7.28),    (7.l8)b),    (7.30)    and    (7.31)   we  have 

(7.331    0<(l-exp)-[    ,         ,?;,   Rp^r',n)   + 

^x   t  Ax    „  )  -1 

for  some  n<N.     However,    by    (7.20a)    and    (7.8) 
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for  some   r"   and  n"<n'<n   .      Thus   K     and  0  are  bounded   for 

all  meshes    A^  satisfying    (5.7)    and    (5.8).      '-^'he    constants 

a,    (3,    y,    5   defined   in   ( 6.  li| ),  >|^  defined  after    (7,1|)    aid 

4   defined   after    (7.3)    are  bounded   independent    of    the  mesh 

spacing.      Thus    the   errors    | e^      |,    which  are    bounded  by 

r,  s 

max  M  of  (7,29),  tend  to  zero  with  the  mesh  spacing, 

n 
^x — >0,  ^u — >0,  6.[i — ^>0  and  ^^2^ — >0;  provided  conditions 

(5,7)  and  (5,8)  are  satisfied.   Hence  the  numerical 

solutions  w'^   converge  uniformly  to  the  solution 
r,  s 

\|/(r, n, u)  in  D^j  on  an  appropriate  sequence  of  meshes  A^, 
The  stability  and  convergence  of  the  proposed  difference 
scheme  are  thus  demonstrated. 
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Figure    I 
Maximum    Energy    Solutions    in    the    Source    Plane 
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Figure    n 
The    Domain    D      with   Indicated    Boundary  and   Initial    Conditions 
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Figure    HI 
Domain   of    Integration    for    i  -  tti    Birtti    Integral     b'Ci/^] 
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